﻿/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Jurgen Van Gael
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
using System;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Decomposition;

namespace dnAnalytics.Statistics.Distributions
{
    /// <summary>
    /// This class implements functionality for the multivariate normal distribution. This distribution
    /// is parameterized by a mean vector and a covariance matrix.
    /// </summary>
    /// <remarks>The distribution will use the <see cref="System.Random"/> by default. 
    /// Users can set the random number generator by using the <see cref="RandomNumberGenerator"/> property.
    /// 
    /// The statistics classes will check all the incoming parameters whether they are in the allowed
    /// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters
    /// to false, all parameter checks can be turned off.</remarks>
    public class VectorNormal
    {
        /// <summary>
        /// The mean of the vector normal distribution.
        /// </summary>
        private readonly Vector mMean;

        /// <summary>
        /// The covariance matrix of the vector normal distribution.
        /// </summary>
        private readonly Matrix mCovariance;

        /// <summary>
        /// The cached Cholesky factorization of the covariance matrix.
        /// </summary>
        private readonly Cholesky mCholesky;

        /// <summary>
        /// Constructs a zero mean vector normal distribution with an identity covariance matrix.
        /// </summary>
        /// <param name="d">The dimension of the vector normal distribution.</param>
        public VectorNormal(int d)
        {
            mMean = new DenseVector(d, 0.0);
            mCovariance = new DenseMatrix(d);
            for (int i = 0; i < d; i++)
            {
                mCovariance[i, i] = 1.0;
            }
            mCholesky = new Cholesky(mCovariance);

            if (Control.CheckDistributionParameters)
            {
                CheckParameters(mMean, mCholesky);
            }

            RandomNumberGenerator = new System.Random();
        }

        /// <summary>
        /// Constructs a vector normal distribution with specific mean and covariance matrices.
        /// </summary>
        /// <param name="mean">The mean of the vector normal.</param>
        /// <param name="covariance">The covariance of the vector normal distribution.</param>
        public VectorNormal(Vector mean, Matrix covariance)
        {
            mMean = mean;
            mCovariance = covariance;
            mCholesky = new Cholesky(mCovariance);

            if (Control.CheckDistributionParameters)
            {
                CheckParameters(mMean, mCholesky);
            }

            RandomNumberGenerator = new System.Random();
        }

        /// <summary>
        /// A string representation of the distribution.
        /// </summary>
        public override string ToString()
        {
            return "VectorNormal(Dimension = " + mMean.Count + ")";
        }

        /// <summary>
        /// Checks the parameters of the vector normal distribution.
        /// </summary>
        /// <param name="mean">The mean of the vector normal.</param>
        /// <param name="covariance">The covariance of the vector normal distribution.</param>
        /// <exception cref="NotPositiveDefiniteException">If the covariance matrix is not positive definite.</exception>
        public static void CheckParameters(Vector mean, Matrix covariance)
        {
            // Try a matrix factorization of S.
            Cholesky chol = new Cholesky(covariance);
            if (!chol.IsPositiveDefinite())
            {
                throw new NotPositiveDefiniteException();
            }
        }

        /// <summary>
        /// Checks the parameters of the vector normal distribution.
        /// </summary>
        /// <param name="mean">The mean of the vector normal.</param>
        /// <param name="chol">The Cholesky factorization of the covariance matrix.</param>
        /// <exception cref="NotPositiveDefiniteException">If the covariance matrix is not positive definite.</exception>
        public static void CheckParameters(Vector mean, Cholesky chol)
        {
            if (!chol.IsPositiveDefinite())
            {
                throw new NotPositiveDefiniteException();
            }
        }

        /// <summary>
        /// The dimensionality of the vector normal distributed variable.
        /// </summary>
        public int Dimension
        {
            get { return mMean.Count; }
        }

        /// <summary>
        /// Gets the mean.
        /// </summary>
        /// <value>The mean.</value>
        public Vector Mean
        {
            get { return mMean; }
        }

        /// <summary>
        /// Gets the covariance.
        /// </summary>
        /// <value>The covariance.</value>
        public Matrix Covariance
        {
            get { return mCovariance; }
        }

        /// <summary>
        /// Gets the mode.
        /// </summary>
        /// <value>The mode.</value>
        public Vector Mode
        {
            get { return mMean; }
        }

        /// <summary>
        /// Gets the median.
        /// </summary>
        /// <value>The median.</value>
        public Vector Median
        {
            get { return mMean; }
        }

        /// <summary>
        /// Gets the entropy.
        /// </summary>
        /// <value>The entropy.</value>
        public double Entropy
        {
            get
            {
                double a = Dimension * Math.Constants.LogSqrt2PiE;
                double b = System.Math.Log(mCholesky.Determinant()) / 2.0;
                return a + b;
            }
        }

        /// <summary>
        /// Evaluates the probability density function for a the vector normal distribution.
        /// </summary>
        public double Density(Vector x)
        {
            Vector d = (x - mMean);
            Vector e = mCholesky.Solve(d);

            return System.Math.Exp(-0.5 * d.DotProduct(e)) / System.Math.Pow(2.0 * Math.Constants.Pi, Dimension / 2.0)
                                                           / System.Math.Sqrt(mCholesky.Determinant());
        }

        /// <summary>
        /// Gets or sets the random number generator.
        /// </summary>
        /// <value>The random number generator used to generate a random sample.</value>
        public System.Random RandomNumberGenerator { get; set; }

        /// <summary>
        /// Samples a vector normal distributed random variable.
        /// </summary>
        public Vector Sample()
        {
            return Sample(RandomNumberGenerator, mMean, mCholesky);
        }

        /// <summary>
        /// Samples a vector normal distributed random variable.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="mean">The mean of the vector normal distribution.</param>
        /// <param name="chol">The Cholesky factorization of the covariance matrix.</param>
        /// <exception cref="NotConformableException">If the number of rows or columns of the covariance matrix are not
        /// the same as the dimension of the mean vector.</exception>
        public static Vector Sample(System.Random rnd, Vector mean, Cholesky chol)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(mean, chol);
            }

            int D = mean.Count;

            // Sample a standard normal variable.
            Vector v = new DenseVector(D, 0.0);
            for (int d = 0; d < D; d += 2)
            {
                double r;
                v[d] = Normal.SampleBoxMuller(rnd, out r);
                if (d + 1 < D)
                {
                    v[d + 1] = r;
                }
            }

            // Return the transformed variable.
            return mean + chol.Factor()*v;
        }

        /// <summary>
        /// Samples a vector normal distributed random variable.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="mean">The mean of the vector normal distribution.</param>
        /// <param name="covariance">The covariance matrix of the vector normal distribution.</param>
        /// <exception cref="NotConformableException">If the number of rows or columns of the covariance matrix are not
        /// the same as the dimension of the mean vector.</exception>
        public static Vector Sample(System.Random rnd, Vector mean, Matrix covariance)
        {
            Cholesky chol = new Cholesky(covariance);
            return Sample(rnd, mean, chol);
        }

        /// <summary>
        /// Samples a set of vector normal distributed random variables.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="n">The number of vectors that need to be drawn.</param>
        /// <param name="mean">The mean of the vector normal distribution.</param>
        /// <param name="chol">The Cholesky factorization of the covariance matrix.</param>
        /// <exception cref="ArgumentOutOfRangeException">If the number of rows or columns of the covariance matrix are not
        /// the same as the dimension of the mean vector.</exception>
        public static Vector[] Sample(System.Random rnd, int n, Vector mean, Cholesky chol)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(mean, chol);
            }

            int D = mean.Count;

            // The Cholesky factor for the covariance matrix.
            Matrix L = chol.Factor();

            // The array of vectors to return.
            Vector[] ret = new Vector[n];
            for (int i = 0; i < n; i++)
            {
                // Sample a standard normal variable.
                Vector v = new DenseVector(D, 0.0);
                for (int d = 0; d < D; d += 2)
                {
                    double r;
                    v[d] = Normal.SampleBoxMuller(rnd, out r);
                    if (d + 1 < n)
                    {
                        v[d + 1] = r;
                    }
                }

                // Transform the variable with mean and covariance.
                ret[i] = mean + L*v;
            }

            return ret;
        }

        /// <summary>
        /// Samples a set of vector normal distributed random variables.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="n">The number of vectors that need to be drawn.</param>
        /// <param name="mean">The mean of the vector normal distribution.</param>
        /// <param name="covariance">The covariance matrix of the vector normal distribution.</param>
        /// <exception cref="ArgumentOutOfRangeException">If the number of rows or columns of the covariance matrix are not
        /// the same as the dimension of the mean vector.</exception>
        public static Vector[] Sample(System.Random rnd, int n, Vector mean, Matrix covariance)
        {
            Cholesky chol = new Cholesky(covariance);
            return Sample(rnd, n, mean, chol);
        }
    }
}